********************************************************************************
* Modelling of representative dataset from Air Products
* Lei Zhang
* 11/11/2013
********************************************************************************

SETS
    i   "Set of Air Products' plants"       /P1*P5/
    ii  "Set of Competitor's plants"        /PC1*PC4/
    j   "Set of Demand points(markets)"     /D1*D20/
    k   "Set of products"                   /O1/
    t   "Set of time periods"               /T1*T67/;

ALIAS(t,tt);

*-------------------------------------------------------------------------------

SCALAR CMAX "Max amount of expansion" /100/;
SCALAR CMIN "Min amount of expansion" /1/;

$include tabledata

* TABLE DE(j,k,t)        "Demand of product each market at each time period(100ton/90day)"
* TABLE DIS(i,j)         "Distance between Air Products' plants and markets(mile)"
* TABLE DISC(ii,j)       "Distance between competitor's plants and markets(mile)"
* TABLE PRICE(k,t)       "Price of products at different time periods(MM$/100ton)"
* TABLE FALPHA(i,k,t)    "Fixed cost(MM$)"
* TABLE FBETA(i,k,t)     "Variable cost of Air Products(MM$/100ton)"
* TABLE FGAMMA(i,k,t)    "Expansion cost of Air Products(MM$/100ton)"
* TABLE FGAMMAH(i,k,t)   "First 100ton higher expansion cost part of Air Products(MM$/100ton)"
* TABLE FT(i,k,t)        "Distribution cost of Air Products(MM$/100ton/mile)"
* TABLE FCALPHA(ii,k,t)  "Fixed cost of competitor(MM$)"
* TABLE FCBETA(ii,k,t)   "Variable cost of competitor(MM$/100ton)"
* TABLE FCGAMMA(ii,k,t)  "Expansion cost of competitor(MM$/100ton)"
* TABLE FCT(ii,k,t)      "Distribution cost of competitor(MM$/100ton/mile)"

FGAMMA(i,k,t) = FGAMMA(i,k,t) / 10;
FGAMMAH(i,k,t) = FGAMMAH(i,k,t) / 10;
FCGAMMA(ii,k,t) = FCGAMMA(ii,k,t) / 10;

*-------------------------------------------------------------------------------

PARAMETER cc(ii);

VARIABLES
    npv             "Net present value"
    cost            "Market cost";

BINARY VARIABLES
    x(i,k,t)        "Selection of capacity Air Products' expansion"
    z(i,k,t)        "If first expansion happend (if c > 0), z=sum(w)"
    w(i,k,t)        "When first expansion happend";

POSITIVE VARIABLES
    y(i,j,k,t)      "Amount of product that Air Products sells to market"
    yc(ii,j,k,t)    "Amount of product that competitor sells to market"
    c(i,k,t)        "Capacity of Air products facility"
    dc(i,k,t)       "Capacity expansion of Air products";

    c.fx('P1','O1','T1') = 72;
    c.l('P1',k,t) = 72;
    c.fx('P2','O1','T1') = 225;
    c.l('P2',k,t) = 225;
    c.fx('P3','O1','T1') = 450;
    c.l('P3',k,t) = 450;
    c.fx('P4','O1','T1') = 0;
    c.l('P4',k,t) = 0;
    c.fx('P5','O1','T1') = 0;
    c.l('P5',k,t) = 0;

    cc('PC1') = 270;
    cc('PC2') = 450;
    cc('PC3') = 90;
    cc('PC4') = 630;

    w.fx('P1',k,'T1') = 1;
    w.fx('P2',k,'T1') = 1;
    w.fx('P3',k,'T1') = 1;

*-------------------------------------------------------------------------------

EQUATIONS
    objout          "Objective function maximize Air Products NPV"
    objin           "Objective function minimize market cost"
    epd(i,k,t)      "Investment decision in Air Products capacity expansion"
    upb(i,k,t)      "Upper bound of capacity expansion"
    lob(i,k,t)      "Lower bound of capacity expansion"
    tmd(j,k,t)      "Demand satisfaction for all markets"
    spd(i,k,t)      "Capacity and supply of Air Products"
    spdc(ii,k,t)    "Capacity and supply of competitor"
    b1(i,k,t)       "When expansion exist, first expansion exist"
    b2(i,k,t)       "relation between z and w"
;

* npv = total_income - fixed_cost - variable_cost - expansion_cost
objout..            npv =e= SUM((i,j,k,t), PRICE(k,t)*y(i,j,k,t))
                        - SUM((i,k,t), FALPHA(i,k,t)*z(i,k,t))
                        - SUM((i,k,t), FBETA(i,k,t)*c(i,k,t))
                        - SUM((i,k,t), FGAMMAH(i,k,t)*w(i,k,t)+FGAMMA(i,k,t)*dc(i,k,t));
* cost = market_cost + transportation_cost;
objin..             cost =e= SUM((i,j,k,t), PRICE(k,t)*y(i,j,k,t))
                        + SUM((ii,j,k,t), PRICE(k,t)*yc(ii,j,k,t))
                        + SUM((i,j,k,t), FT(i,k,t)*DIS(i,j)*y(i,j,k,t))
                        + SUM((ii,j,k,t), FCT(ii,k,t)*DISC(ii,j)*yc(ii,j,k,t));

* c(i,k,t) = c(i,k,t-1) + dc(i,k,t) (when t < card(t)-1)
epd(i,k,t)$(ord(t)+1 le card(t)).. c(i,k,t+1) =e= c(i,k,t) + dc(i,k,t);
* dc(i,k,t) <= CMAX*x(i,k,t)
upb(i,k,t)..        dc(i,k,t) =l= CMAX*x(i,k,t);
* dc(i,k,t) >= CMIN*x(i,k,t)
lob(i,k,t)..        dc(i,k,t) =g= CMIN*x(i,k,t);
* sumi_y + sumi_yc = market_j's_demand
tmd(j,k,t)..        SUM(i, y(i,j,k,t)) + SUM(ii, yc(ii,j,k,t)) =e= DE(j,k,t);
* sumj_y < c(i,k,t)
spd(i,k,t)..        SUM(j, y(i,j,k,t)) =l= c(i,k,t);
* sumj_yc < cc(ii,k,t)
spdc(ii,k,t)..      SUM(j, yc(ii,j,k,t)) =l= cc(ii);

b1(i,k,t)..         x(i,k,t) =l= z(i,k,t);
b2(i,k,t)..         z(i,k,t) =e= SUM(tt$(ord(t) ge ord(tt)), w(i,k,tt));

*-------------------------------------------------------------------------------

MODEL COMPETITION /all/;

$echo bilevel y c dc x z w min cost objin tmd spdc > "%emp.info%";

SOLVE COMPETITION USING emp MAXIMIZING npv;
